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We present a calculation of a fourth-order, time-dependent density correlation function that mea- 
sures higher-order spatiotemporall correlations of the density of a liquid. From molecular dynamics 
simulations of a glass-forming Lennard- Jones liquid, we find that the characteristic length scale of 
this function has a maximum as a function of time which increases steadily beyond the characteristic 
length of the static pair correlation function g{r) in the temperature range approaching the mode 
coupling temperature from above. 

PACS numbers: PACS numbers: 64.70.pf, 61.20.1c 



Relaxation in liquids near their glass transition involves the correlated motion of groups of neighboring particles j^, 
^, ^. This correlated motion results in spatially heterogeneous dynamics, which becomes increasingly heterogeneous 
as the liquid is cooled. Much remains to be understood regarding the nature of this heterogeneity, and how to best 
measure and quantify it. The traditional two-point, time-dependent, van Hove density correlation function G{r,t), 
provides information about the transient "caging" of particles on cooling |Q] , but does not provide local information 
about correlated motion and dynamical heterogeneity. In particular, the static correlation length associated with two- 
point density fluctuations remains relatively constant upon cooling j^]. Instead, other correlation functions, which 
involve, e.g., spatial correlations of the displacements of particles in the liquid, and other measures of correlated 
motion, have been used to demonstrate the heterogeneous nature of the liquid dynamics in molecular dynamics 
(MD) simulations |Q, ^. These "measures" are readily accessible in colloidal suspensions, where microscopy provides 
detailed information on particle trajectories similar to the information obtained from MD simulations. Indeed, such 
experiments have confirmed simulation predictions of increasingly heterogeneous dynamics near the glass transition 

In this paper, we evaluate a fourth-order, time-dependent density correlation function g4(r, t) that is more sensitive 
to spatially heterogeneous dynamics than G{r, t). This four-point function was first investigated in supercooled liquids 
by Dasgupta et al. |l0[] , but they did not detect a growing correlation length in their simulations. Recently, it was 
shown ||ll, |l2j that it is possible to define a generalized, time-dependent susceptibility X4(^) which (i) is proportional 
to the volume integral of 54(7", t) in the same way that the isothermal compressibility is related to the volume integral 
of the static pair correlation function [ p^ , (ii) is non-zero in the caging regime, and (iii) increases with decreasing 
T. While this indirectly suggests the presence of a growing correlation length, neither g4{r,t) nor its range ^4{t) was 
calculated in those works. Here we calculate gi{r,t) in a simulation of 8000 Lennard-Jones (LJ) particles, and show 
that £,4:{t) grows slowly but steadily beyond the correlation length ^ of the static pair correlation function g{r) as 
temperature T is decreased from the onset of caging towards the mode coupling temperature Tmct 0- 

We first briefly review the general theoretical framework, some of which was previously discussed in Refs. jlll, |lj, 
and extend it to obtain a form for gi{r, t) suitable for calculation in our simulations. Consider a liquid of N particles 
occupying a volume V , with density p(r, t) = X^iLi ^{'^ ~ Extending an idea originally proposed for spin glasses 

one may construct a time-dependent 'order parameter' that compares the liquid configuration at two different 
times IT] 
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Q{t)= /'dridr2p(ri,0)p(r2,t)u;(|ri-r2|)=^^«;(Ir,(0)-r,(t)|). (1) 
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Here, in the second equality refers to the position of particle i, and it;(|ri — rj]) is an "overlap" function which is 
unity for jr^ — rJ < a and zero otherwise, where the parameter a is associated with the typical vibrational amplitude 
of the particles |ll|, For the present work, we take a — 0.3 particle diameters as in Refs. [HJ 0, since this value 
maximizes the effect studied here As defined, Q{t) is the number of "overlapping" particles when configurations 
of the system at i = and at a later time t are compared; that is, Q(t) counts the number of particles that cither 
remain within a distance a of their original position, or are replaced (within a distance a) by another particle in an 
interval t. 

The fluctuations in Q{t) are described by a generalized susceptibility 

Xiit) = ^{{Q'{t))^{Qit)y), (2) 

where /3 = l/fcsT, ks is Boltzmann's constant, and (■ • •) indicates an ensemble average as in Ref. Note that 

at very early times, when Q{t) — N because no particle has yet moved beyond a distance a, Xiit) is identical to the 
isothermal compressibility kt 10 • Substituting Eq. (|l|) into Eq. (|2|), we obtain 



Xiit) 

where 
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j dvidv2dv3dri G4(ri, r2, rs, r4, i), (3) 



G'4(ri,r2,r3,r4,t) = (p(ri, 0)p(r2, i)M;(|ri - r2|)p(r3, 0)p(r4, t)w(|r3 - r4|)) 

- (p(ri,0)p(r2,tM|ri-r2|))(p(r3,0)p(r4,i)i«(|r3-r4|)). (4) 

As written, G/^ is a function of four spatial and two temporal coordinates, constrained in a particular way by the 
overlap functions. To investigate the behavior of 6*4, we consider a function g^^r^t) such that X4{t) = P J dr gi(r,t). 
We integrate first over V2 and r4 in Eq. (||) , define r = r3 — ri , and then integrate over r3 to obtain 

54(r-, t)^{^ yi - ^fc(O) + r.(0))u;(|r,(0) - v,{t)\)<\^m - - (^)' (5) 

^ ijkl 

We investigate the behavior of g4{r,t), which is the angular averaged function of a single variable r. With the above 
choice of integration variables, gi{r,t) describes spatial correlations between overlapping particles at the initial time 
(using information at time t to label the overlapping particles). The first term in (?4(r, t) is a pair correlation function 
restricted to the subset of overlapping particles, gf{r,t). The second term represents the "bulk" probability of any 
two particles overlapping. We can rewrite g^i^r^t) as 



,oVwA /Q{t)\^ _/Q{t)\^\gf{r,t) 



g,ir,t) ^ gfir,t)-{^:^) = (^:^) 
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- {^)\lir.t). (6) 

Since (^(p) is a function of time, information about spatial correlations is contained in gX{r, t), which we investigate 
in the rest of the paper. In comparing 34 (r, t) with other correlation functions used to study glass-forming liquids, we 
note that neutron scattering studies, such as those by Colmenero and Richter and coworkers on polymer systems [ p^ , 
measure at most two-point spatiotemporal density correlation functions. 4-D NMR methods used to probe dynamical 
heterogeneity measure multiple time correlation functions p9| . In contrast to those, 54 (r, t) contains additional, higher 
order spatiotemporal information. 

To evaluate gl{r,t), we perform MD simulations of a model LJ glass- forming liquid. The system is a three- 
dimensional (50:50) binary mixture of 8000 particles interacting via LJ interaction parameters |2^. We simulate 
eight state points in the microcanonical ensemble at fixed density p = 1.29, in a temperature range 0.59 — 2.0; we 
estimate Tmct = 0.57 ± 0.02 and the Vogel-Fulcher temperature Tg = 0.48 ± 0.02. The error bars are confidence 
intervals obtained as a result of fitting Tq, (T) to a power law and exponential form, respectively. The onset of caging 
and non-exponential relaxation of the intermediate scattering function occur at Tcago ~ 1-0. Our simulations thus span 
a range that includes the onset of caging and the initial slowing down of the liquid approaching Tmct. The isochore 
along which the simulations are performed was chosen to reproduce simulations of a smaller system size performed 
earlier, which showed (i) that this model exhibits spatially heterogeneous dynamics||l2|, and (ii) that transitions 
between inherent structures close to Tmct occur through the collective, quasi-one-dimensional motion of strings of 
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particles | |2l[ |. All data evaluated in the present work are in thermodynamic equilibrium above the glass transition 
temperature. 

We plot gl{r, t) for several t at our second coldest temperature T = 0.60 in Figs. |l|(a) and (b). The positions of the 
peaks in g\{r,t) are identical to the positions of the peaks in gir) (not shown). We confirm that gl{r,t) — g{r) — 1 
in the ballistic and diffusive regime. Note that in the diffusive regime g°^{r,t) is the pair correlation function of the 
random overlaps normalized by (^^)^ to yield g{r). We plot Xiit) and ( ^ffi ) at T = 0.60 in Fig. 0(c); as found in 
Refs. p] , [T^ , Xiit) has a maximum at an intermediate time t*l. Rcfs. 12 showed that this maximum increases 
and shifts to longer time as T decreases toward Tmct , and we further find that is in the a-relaxation regime at each 
T. gl{r, t) deviates from g{r) — 1 when (^^) deviates from unity and Xii^) becomes non-zero|^^. The amplitude and 
range of g|(r, t) increase with increasing t until a time t\. At t\, gl{r, t|) (indicated by the solid curve in Figs. |^(a,b)) 
exhibits a long tail which decreases slowly to zero with increasing distance. For t greater than tj, the amplitude and 
range of gl{r,t) decrease, and gl{r,t) = g{r) — 1 when Xii^) decays to zero (not shown). The positions of the peaks 
in gl{rA) do not appear to change with decreasing T. 

Fig. H shows the T dependence of gUr.t) at the peak characteristic time tl when the correlations at each T are 
most pronounced, as measured by Xiit)- The inset of Fig. || shows the "four-point" structure factor Sl{q,tl) ~ 
p J drgl{r,tl)sin{qr)/qr, and static structure factor S{q) — l=pj dr[g{r) — 1] sin(gr)/qr. We find that while S{q) 
shows no change at small g, Si{q,t^) develops a peak at small q that grows with decreasing T, indicating a growing 
range of correlations. 

Fig. ^ provides a close-up of the behavior of gl{r, tl) for 1.7 < r < 7, for several values of T. To extract a value for 
^4(4), we fit peaks of gl{r,t) in the range shown to the exponential function y{r) = a * exp(— r/^4(t)). We refer to this 
method as an "envelope fit". The time dependence of ^4{t) obtained from this fit is plotted for several state points in 
Fig. 1^. We see that the qualitative behavior of ^4(1) is similar to that of Xiit)- ^ii^) has a maximum in time, and as 
T decreases, the amplitude and time of this maximum increase. 

The length scale ^4(t) characterizes the typical distance over which "overlapping" particles are spatially correlated, 
and includes contributions from the static two-point density correlations. At temperatures above the onset of caging, 
where the dynamics is everywhere homogeneous, ^4(^4) and ^ coincide. Below the onset of caging, ^4(^4) begins to 
grow larger than ^; over the limited T-range of our simulations, we find that ^4{t) increases from 0.8 ± 0.1 particle 
diameters above Tcage to 1.5 ± 0.1 particle diameters within 5% of Tmct (Fig- § inset). In contrast, we find that 
^, which is coincident with $,4{t) at short times when (^j^) = 1, changes less, from 0.8 ± 0.1 to 1.1 ± 0.1 particle 
diameters. Thus, over the T-range studied, the characteristic distance over which particle motion is most correlated 
grows to exceed the static correlation length. While we cannot reliably predict the behavior of S,4(tl) at lower T, we 
find no tendency for slowing down of its growth, unlike the low-T behavior of ^ which is known to remain finite and 
small. We do observe a slight depression of ^4{t) and X4(i) at our lowest T = 0.59 (not shown), but we speculate this 
may be due to finite size effects, which we will explore in detail in a later publication. 

The relatively small but growing correlation length calculated here from the four-point spatiotemporal density 
correlation function should be contrasted with that characterizing the size of highly mobile regions within the fluid 
0. That length was shown to grow much more rapidly on cooling, approaching the size of the simulation box close 
to Tmct- Interestingly, whereas the correlation length of highly mobile regions was found to be largest on a time 
scale in the /^-relaxation regime, the length calculated in the present paper is largest in the a-regime. We note that 
g4{r,t) is dominated by "caged" particles, and thus ^4(t) may be related to length scales calculated in Rcf. The 
relationship of the different length scales characterizing dynamical heterogeneity will be explored in a subsequent 
publication. 
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FIG. 1: Time dependence of gl{r,t) at T = 0.60, showing (a) the amplitude and range of gl{r,t) growing in time and (b) the 
amphtude and range of gl{r,t) decaying in time. The fraction (■^^) indicates the average fraction of overlapping particles 
present at time t. (c) Time dependence of X4(^) a-nd (^j^) at T = 0.60. The long time value of (^^) = pVa (solid line), 

where Va = ^^y^ corresponds to the probability of finding a random overlapping particle. The error bars in Q(t) and Xi{t) are 
calculated by averaging over 100 successive independent blocks, and are representative of the error in all points. 




FIG. 2: Temperature dependence of gl{r,t) at six values of T, indicated in legend. Insets show the structure factor S^q, 
and static structure factor S{q) for the same values of T. We note that g2 (r, t) is analogous to p(r) — 1, and SI {q, t) is analog 
to S{q) - 1. 
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FIG. 3: Somilogarithmic plot of 54 (r, t^) in the range 1.7 < r < 7 at four T. The soUd hnes are the exponential curves obtained 
from an "envelope fit". 




FIG. 4: Time and T dependence of (,A{t). The inset shows iA{ti) vs. T. The data shown is smoothed by performing running 
average over successive groups of five points. 



